The Effect of Blood Rheology and Inlet Boundary Conditions on Realistic Abdominal Aortic Aneurysms under Pulsatile Flow Conditions

Background: The effects of non-Newtonian rheology and boundary conditions on various pathophysiologies have been studied quite extensively in the literature. The majority of results present qualitative and/or quantitative conclusions that are not thoroughly assessed from a statistical perspective. Methods: The finite volume method was employed for the numerical simulation of seven patient-specific abdominal aortic aneurysms. For each case, five rheological models and three inlet velocity boundary conditions were considered. Outlier- and heteroscedasticity-robust ANOVA tests assessed the simultaneous effect of rheological specifications and boundary conditions on fourteen variables that capture important characteristics of vascular flows. Results: The selection of inlet velocity profiles appears as a more critical factor relative to rheological specifications, especially regarding differences in the oscillatory characteristics of computed flows. Response variables that relate to the average tangential force on the wall over the entire cycle do not differ significantly across alternative factor levels, as long as one focuses on non-Newtonian specifications. Conclusions: The two factors, namely blood rheological models and inlet velocity boundary condition, exert additive effects on variables that characterize vascular flows, with negligible interaction effects. Regarding thrombus-prone conditions, the Plug inlet profile offers an advantageous hemodynamic configuration with respect to the other two profiles.


Introduction
Cardiovascular diseases (CVDs) are the leading cause of death globally [1]. An estimated 17.9 million individuals died from CVDs in 2019, representing more than 30% of all deaths worldwide. They comprise a group of disorders that affect the heart and blood vessels such as rheumatic and congenital heart disease, cerebrovascular disease (stroke), aortic aneurysms, and pulmonary embolism. Specifically, abdominal aortic aneurysms (AAAs) represent the 13th leading cause of death in Western societies [2]. It is typically defined as an at least 1.5-fold increase in vessel diameter and its rupture constitutes the most dangerous implication, which is usually accompanied by a catastrophic insult with an overall mortality between 70% and 90% [3,4]. The traditional criterion for rupture risk assessment is the AAA's maximum diameter that should not exceed the threshold of 55 mm. Even though this criterion has served the medical community over the years, it is reasonable to question the concept that a single parameter can sufficiently fulfill the needs of all patients. To this end, additional procedures have been developed towards increasing rupture risk assessment, such as surgery outcome prediction [5,6], noninvasive imaging [7,8], and numerical simulations through computational fluid dynamics (CFD) [9][10][11][12].
The latter approach provides an advantage to researchers and clinicians, since it can resolve spatial and temporal resolution of blood flow in high detail without the need of in vivo measurements, which are invasive and possibly demanding to perform. A prerequisite, though, of almost any CFD simulation is an appropriate set of conditions that formulate a well-defined problem to be solved. These conditions depend on the question under consideration, but typically involve assumptions on blood flow rheology, inlet/outlet boundary or initial conditions, and biomechanical properties of relevant vessels amongst others. To this end, numerous works have appeared in the literature studying the effects of various assumptions on the hemodynamic behavior of idealized or realistic geometries.
Neofytou et al. [13] studied the effects of different blood rheological models on a stenosis and an abdominal aortic aneurysm, focusing on the distribution of wall shear stress in the vicinity of vessel abnormalities. Arzani [14] proposed a residence time non-Newtonian model that accounts for the rouleaux formation time-scale in blood shear-thinning behavior, while Bilgi et al. [15] concluded that a Carreau fluid through a hyperelastic vessel behaves substantially different than a Newtonian fluid with a linearly elastic arterial wall. Skiadopoulos et al. [16] compared three blood flow models in patient-specific cardiovascular systems and concluded that the Newtonian assumption is valid only for high shear and flow rates.
The impact of boundary conditions has also been examined quite extensively. Morbiducci et al. [17] studied idealized versus measured velocity profiles as inlet boundary conditions in the human aorta. Similarly, Youssefi et al. [18] performed a numerical study on the effect of various inlet boundary conditions in the thoracic aorta and concluded that idealized velocity profiles can potentially lead to significant alterations of velocity patterns and magnitudes in the aorta. More studies on aortic flows include the works of Madhavan et al. [19] and Fuchs et al. [20], where it was shown that simulation results were in general sensitive to the choice of boundary conditions. The impact of inlet boundary conditions on blood flow has also been examined in various other physiologies such as stented coronary arteries [21], carotid bifurcations [22], intracranial aneurysms [23], and abdominal aortic aneurysms [24].
The plethora of available studies in the literature present a qualitative and/or quantitative evaluation of the effects of rheological models and boundary conditions on the hemodynamic behavior of various pathophysiologies, without usually assessing observed differences from a statistical perspective. The aim of this work is to analyze the effects of blood flow models and inlet boundary conditions on patient-specific abdominal aortic aneurysms, with statistical tools which are robust to deviations from the standard assumptions of conventional analysis of variance (ANOVA): equal levels of uncertainty across groups (homoscedasticity) and absence of outliers. Specifically, seven AAA cases were considered, with five different models and three inlet boundary conditions, yielding 105 numerical simulations in total. Fourteen response variables that characterize vascular flows were computed per simulation and analyzed with modern statistical methods. The response variables quantify (a) the average tangential force on the wall over the entire cycle; (b) the oscillatory nature of the flow; (c) flow asymmetry; (d) flow dispersion; and (e) the extent of thrombogenic stimulating environments.
The manuscript is organized as follows: Section 2 reports the segmentation/reconstruction and meshing process for the seven cases considered. It also formulates the mathematical framework that underlies numerical simulations and presents the adopted velocity inlet boundary conditions and the alternative rheological specifications. Furthermore, it defines the hemodynamic variables, as well as the robust statistical models and hypothesis tests employed for their analysis. Section 3 presents a mesh convergence study and reports the outcomes of the statistical analysis, which evaluates the effect of inlet velocity and blood flow model on the overall hemodynamic behavior of the AAAs under consideration. Finally, Section 4 presents a thorough discussion on our findings and Section 5 concludes the paper.

Image Segmentation, Surface Reconstruction, and Meshing
Seven AAA patients (denoted 2B, 7A, 14B, 16A, 31A, 41B, and 63A) underwent electrocardiogram-gated (ECG-gated) computed tomography scans to image the AAA geometry. Following [25], a Somatom definition flash, dual-source-dual-energy CT scanner (Siemens, Munich, Germany), with before and after contrast media administration with retrospective ECG-gated spiral acquisition, was used for imaging. A total effective dose of 5.5 mSv at 80 bpm non-ionic contrast was used with a slice thickness of 0.625 mm and image matrix size 512 × 512. The temporal resolution was set equal to 83 ms with 0.33 mm in-plane spatial resolution. Segmentation and surface reconstruction were performed manually using the open source software ITK_SNAP [26]. The generated surfaces were smoothed using VMTK [27]; specifically, the Taubin algorithm that preserves the volume under consideration [28] was implemented. A passband of 0.01 and 100 iterations were sufficient in all cases in order to remove the surface noise, without changing the AAA sac volume more than 0.15%. Finally, cylindrical extensions were added at the aorta and iliac arteries to allow flow development.
All surfaces were then meshed with ANSA (BETA CAE Systems S.A.) using a pure hexahedral mesh. As shown by De Santis et al. [29], hexahedral meshes should be preferred to other types of meshes since they require a fewer number of elements for a specific level of accuracy. An appropriate O-Grid (with 0.25 parametric and 0.95 Bell shape values) was constructed in all cases to capture the high velocity gradients due to the boundary layer close to the rigid wall. Figure 1 presents the inlet (A) and part of the surface (B) meshes.
consideration. Finally, Section 4 presents a thorough discussion on our findings and Sec tion 5 concludes the paper.

Image Segmentation, Surface Reconstruction, and Meshing
Seven AAA patients (denoted 2B, 7A, 14B, 16A, 31A, 41B, and 63A) underwent elec trocardiogram-gated (ECG-gated) computed tomography scans to image the AAA geom etry. Following [25], a Somatom definition flash, dual-source-dual-energy CT scanne (Siemens, Munich, Germany), with before and after contrast media administration wit retrospective ECG-gated spiral acquisition, was used for imaging. A total effective dos of 5.5 mSv at 80 bpm non-ionic contrast was used with a slice thickness of 0.625 mm an image matrix size 512 × 512. The temporal resolution was set equal to 83 ms with 0.33 mm in-plane spatial resolution. Segmentation and surface reconstruction were performed manually using the open source software ITK_SNAP [26]. The generated surfaces wer smoothed using VMTK [27]; specifically, the Taubin algorithm that preserves the volum under consideration [28] was implemented. A passband of 0.01 and 100 iterations wer sufficient in all cases in order to remove the surface noise, without changing the AAA sa volume more than 0.15%. Finally, cylindrical extensions were added at the aorta and ilia arteries to allow flow development.
All surfaces were then meshed with ANSA (BETA CAE Systems S.A.) using a pur hexahedral mesh. As shown by De Santis et al. [29], hexahedral meshes should be pre ferred to other types of meshes since they require a fewer number of elements for a specifi level of accuracy. An appropriate O-Grid (with 0.25 parametric and 0.95 Bell shape values was constructed in all cases to capture the high velocity gradients due to the boundar layer close to the rigid wall. Figure 1 presents the inlet (A) and part of the surface (B meshes.

Simulation Setup, Boundary Conditions, and Rheology Models
Blood was modeled as an incompressible, homogeneous, and non-Newtonian fluid The fluid domain was governed by the coupled system of Navier-Stokes and continuit equations formulated as ( ) where the velocity and pressure fields , P U , respectively, depended on blood density and stress tensor τ . For the needs of this study, the stress tensor was expressed in term of the rate-of-deformation tensor D and the shear rate γ  as follows

Simulation Setup, Boundary Conditions, and Rheology Models
Blood was modeled as an incompressible, homogeneous, and non-Newtonian fluid. The fluid domain was governed by the coupled system of Navier-Stokes and continuity equations formulated as where the velocity and pressure fields U,P, respectively, depended on blood density ρ and stress tensor τ. For the needs of this study, the stress tensor was expressed in terms of the rate-of-deformation tensor D and the shear rate . γ as follows All vascular flows were simulated using a commercial finite volume solver (Fluent 17.2, ANSYS Inc., Canonsburg, PA, USA) with a default criterion for convergence equal to 10 −4 . The SIMPLEC algorithm was chosen for pressure-velocity coupling and a fixed time step of 0.005 s was adopted for a cardiac cycle of 1 s. In order to ensure that all transient effects were washed-out, four cycles were simulated before results were collected. All simulations assumed a rigid wall and the no-slip condition was prescribed at the wall boundary. Furthermore, a transient inlet velocity and outlet pressure were prescribed in all simulations; both profiles ( Figure 2) closely follow Olufsen et al. [30].
Bioengineering 2023, 10, x FOR PEER REVIEW 4 of All vascular flows were simulated using a commercial finite volume solver (Flue 17.2, ANSYS Inc., Canonsburg, (PA), USA) with a default criterion for convergence equ to 10 −4 . The SIMPLEC algorithm was chosen for pressure-velocity coupling and a fix time step of 0.005 s was adopted for a cardiac cycle of 1 s. In order to ensure that all tra sient effects were washed-out, four cycles were simulated before results were collecte All simulations assumed a rigid wall and the no-slip condition was prescribed at the w boundary. Furthermore, a transient inlet velocity and outlet pressure were prescribed all simulations; both profiles ( Figure 2) closely follow Olufsen et al. [30]. The pressure distribution was kept constant for each time step at the outle Throughout the cardiac cycle, though, the distribution followed the one shown in Figu 2b. Regarding velocity, three different profiles (Figure 3), namely Parabolic (Par), Pl (Plug), and Womersley (Wom) were prescribed at the inlet, with equal mean values for three of them ( Figure 2a). This allows investigation of the effect of inlet velocity distrib tion at the hemodynamic behavior of patient-specific aneurysm geometries, which cons tutes the first goal of this work. A user-defined function (UDF) was implemented into the solver to assign the prop velocity according to the following expressions. The pressure distribution was kept constant for each time step at the outlets. Throughout the cardiac cycle, though, the distribution followed the one shown in Figure 2b. Regarding velocity, three different profiles (Figure 3), namely Parabolic (Par), Plug (Plug), and Womersley (Wom) were prescribed at the inlet, with equal mean values for all three of them ( Figure 2a). This allows investigation of the effect of inlet velocity distribution at the hemodynamic behavior of patient-specific aneurysm geometries, which constitutes the first goal of this work. All vascular flows were simulated using a commercial finite volume solver (Fluent 17.2, ANSYS Inc., Canonsburg, (PA), USA) with a default criterion for convergence equal to 10 −4 . The SIMPLEC algorithm was chosen for pressure-velocity coupling and a fixed time step of 0.005 s was adopted for a cardiac cycle of 1 s. In order to ensure that all transient effects were washed-out, four cycles were simulated before results were collected. All simulations assumed a rigid wall and the no-slip condition was prescribed at the wall boundary. Furthermore, a transient inlet velocity and outlet pressure were prescribed in all simulations; both profiles ( Figure 2) closely follow Olufsen et al. [30]. The pressure distribution was kept constant for each time step at the outlets. Throughout the cardiac cycle, though, the distribution followed the one shown in Figure  2b. Regarding velocity, three different profiles (Figure 3), namely Parabolic (Par), Plug (Plug), and Womersley (Wom) were prescribed at the inlet, with equal mean values for all three of them ( Figure 2a). This allows investigation of the effect of inlet velocity distribution at the hemodynamic behavior of patient-specific aneurysm geometries, which constitutes the first goal of this work. A user-defined function (UDF) was implemented into the solver to assign the proper velocity according to the following expressions. A user-defined function (UDF) was implemented into the solver to assign the proper velocity according to the following expressions.
Equation (3) defines the streamwise velocity distributions at the inlet in terms of the radius, R, angular frequency, ω, pressure gradient along the flow, ∂P/∂z, number of Fourier modes, n, Womersley parameter, α = R √ ωρ/µ, and Bessel function, J 0 . For the needs of this study, fourteen modes were used in the Fourier series expansion to ensure that deviations of the constructed profiles did not deviate more than 0.5% with respect to the ones presented in [30]. In all cases, the flow is assumed laminar with a mean and peak Reynolds number as well as Womersley parameter, α, presented in Table 1. This assumption is in accordance with [31][32][33][34][35], where it was assumed that blood flow is laminar, even during exercise, for asymmetric AAAs [36]. Care must be taken, though, since turbulence might emerge under certain circumstances. As pointed out in [32], severe angulation of the proximal neck can cause strong turbulence, altering the distribution of WSS on the artery wall. Additionally, Khanafer et al. [37] showed that the increased shear stress due to local turbulence can generate further dilation of the aneurysm sac, creating a mechanism for aneurysmal growth and potential rupture.
The second research question that is investigated herein evaluates the effect of the adopted rheological model on vascular flows. Five specifications were implemented ( Table 2); it should be stressed that for shear rate values greater than 100 s −1 , all models converge to the Newtonian case. Specifications were selected following a recent classification of 16 rheology models [38], which revealed a partition in three main homogeneous groups (clusters) and six in total, with the Newtonian model appearing as an outlier. The Carreau-Yasuda and the Casson models were members of the largest cluster, whereas the Power law was the best representative of the second largest cluster and Herschel-Bulkley a satisfactory representative of the third cluster. Table 2. Rheology models and parameter values. [14,39,40] Casson (Cs) µ [43,44] Power law (P) µ

Hemodynamic and Flow Parameters
Various parameters have been used in the literature to study vascular flows in terms of the wall shear stress vector (WWS) over the cardiac cycle, T. These include the time average wall shear stress (TAWSS, Pa), oscillatory shear index (OSI), and relative residence time (RRT, Pa −1 ), defined by the following expressions According to the above definitions, TAWSS quantifies the average tangential force on the wall over the entire cycle but does not provide any insight on the oscillatory nature of the flow. To this end, He et al. [48] introduced the non-dimensional parameter OSI in order to capture the cyclic variation in WSS. As a result, uniaxial flows yield OSI = 0 while flows with no preferred direction correspond to OSI = 0.5. Finally, Himburg et al. [49] presented RRT to quantify the time that blood resides close to the wall while accounting for the effect of both TAWSS and OSI. The abovementioned variables provide valuable information for various diseased states, such as thrombogenic stimulating environments for TAWSS < 0.4 Pa [50], OSI > 0.3 [49], and RRT > 10 Pa −1 [51].
An additional set of parameters that characterize flows are flow asymmetry, f A , and flow dispersion, f D ; both quantities were calculated on random planes (to ensure that results were not affected in a systematic way) of the aneurysm sac and quantify the eccentricity and broadness of flows, respectively. Following [18,52], the centroid coordinates x 0 , y 0 , z 0 of the top 15% peak systolic velocity V 15% max were compared with respect to the centroid x c , y c , z c of the plane under consideration. Flow asymmetry is formulated as where R eq is the equivalent radius of a circle with the same area as the plane vessel. Thus, flows with low asymmetry values do not deviate significantly from the plane centroid, while flows with high values are eccentric and develop mainly close to the vessel wall. In a similar manner, yielding a broad velocity profile for large dispersion values and a sharper one as values decrease. Figure 4 depicts the seven patient-specific geometries and the corresponding planes where flow asymmetry and dispersion were calculated.
Bioengineering 2023, 10, x FOR PEER REVIEW 7 of 21 yielding a broad velocity profile for large dispersion values and a sharper one as values decrease. Figure 4 depicts the seven patient-specific geometries and the corresponding planes where flow asymmetry and dispersion were calculated.

Statistical Analysis
The effects of alternative inlet velocity distributions ( Figure 3) and rheological models (Table 2) at the hemodynamic behavior of real aneurysm geometries was investigated with two-way ANOVA tests, which evaluate simultaneously the effect of two grouping variables (or factors) on a response variable [53]. The level combinations of the two grouping variables, denoted by IVD and Rm, respectively, produced a balanced design, as the sample sizes within each level of the independent factors are equal. Fourteen response variables that characterize vascular flows were examined: flow asymmetry (FApct), flow dispersion (FDpct), percentage areas with thrombus-prone conditions, designated by the levels of TAWSS (TAWSSpct: % area with TAWSS < 0.4 Pa), OSI (OSIpct: % area with OSI > 0.3), and RRT (RRTpct: % area with RRT > 10 Pa −1 ), and the observed average levels, minima and maxima of hemodynamic variables, namely TAWSS (TAWSSave, TAWSSmin, TAWSSmax), OSI (OSIave, OSImin, OSImax), and RRT (RRTave, RRTmin, RRTmax).
Two-way ANOVA designs evaluate the following initial (null) hypotheses per response variable [53]: (a) there are no significant differences in the observed means across the three alternative IVD choices; (b) there are no significant differences in the observed means across the five alternative rheological models; (c) there is no significant interaction between IVD and Rm. If an ANOVA test signifies substantial evidence against the null hypothesis of equal means across factor levels, the analysis proceeds in the second stage, which comprises multiple pairwise comparisons between the group means: the latter determine if specific group pairs (for Rm, say P versus N, or for IVD, Plug versus Parabolic) are significantly different. On the other hand, if the p-values from an ANOVA test do not provide evidence against the null hypothesis of equal means across groups, there is no need to conduct a post hoc test to determine which groups are different from each other.
Conventional ANOVA is based on group means and assumes normality and homoscedasticity (equal variances across groups) for model residuals. Violation of the abovementioned assumptions may yield inaccurate confidence intervals and poor characterization regarding the significance of observed group differences [54]. Preliminary analyses of the data presented in the next section revealed strong evidence against the assumption of homoscedastic residuals. To safeguard quantitative analyses from false hypotheses, we adopted trimmed-mean-based methods, which are robust to deviations from normality and homoscedasticity. Specifically, the robust two-way ANOVA procedure and the corresponding post hoc tests presented in (p. 335, [54]) were implemented, using R package

Statistical Analysis
The effects of alternative inlet velocity distributions ( Figure 3) and rheological models ( Table 2) at the hemodynamic behavior of real aneurysm geometries was investigated with two-way ANOVA tests, which evaluate simultaneously the effect of two grouping variables (or factors) on a response variable [53]. The level combinations of the two grouping variables, denoted by IVD and R m , respectively, produced a balanced design, as the sample sizes within each level of the independent factors are equal. Fourteen response variables that characterize vascular flows were examined: flow asymmetry (FA pct ), flow dispersion (FD pct ), percentage areas with thrombus-prone conditions, designated by the levels of TAWSS (TAWSS pct : % area with TAWSS < 0.4 Pa), OSI (OSI pct : % area with OSI > 0.3), and RRT (RRT pct : % area with RRT > 10 Pa −1 ), and the observed average levels, minima and maxima of hemodynamic variables, namely TAWSS (TAWSS ave , TAWSS min , TAWSS max ), OSI (OSI ave , OSI min , OSI max ), and RRT (RRT ave , RRT min , RRT max ).
Two-way ANOVA designs evaluate the following initial (null) hypotheses per response variable [53]: (a) there are no significant differences in the observed means across the three alternative IVD choices; (b) there are no significant differences in the observed means across the five alternative rheological models; (c) there is no significant interaction between IVD and R m . If an ANOVA test signifies substantial evidence against the null hypothesis of equal means across factor levels, the analysis proceeds in the second stage, which comprises multiple pairwise comparisons between the group means: the latter determine if specific group pairs (for R m , say P versus N, or for IVD, Plug versus Parabolic) are significantly different. On the other hand, if the p-values from an ANOVA test do not provide evidence against the null hypothesis of equal means across groups, there is no need to conduct a post hoc test to determine which groups are different from each other.
Conventional ANOVA is based on group means and assumes normality and homoscedasticity (equal variances across groups) for model residuals. Violation of the abovementioned assumptions may yield inaccurate confidence intervals and poor characterization regarding the significance of observed group differences [54]. Preliminary analyses of the data presented in the next section revealed strong evidence against the assumption of homoscedastic residuals. To safeguard quantitative analyses from false hypotheses, we adopted trimmed-mean-based methods, which are robust to deviations from normality and homoscedasticity. Specifically, the robust two-way ANOVA procedure and the corresponding post hoc tests presented in (p. 335, [54]) were implemented, using R package WRS2 [55]. In the post hoc comparisons, both confidence intervals and p-values were adjusted to account for multiple testing by controlling the family-wise error rate (or the probability of false significance), using the Hochberg method [54].
The null hypotheses in the statistical tests imply that the main characteristics of vascular flows do not differ substantially for alternative IVD and R m choices. In what follows, to evaluate evidence that favors the alternative hypothesis, we use p-values following the recent recommendation in [56]: values between 0.005 and 0.05 offer weak or "suggestive" evidence, whereas values lower than 0.005 provide strong evidence against the null hypothesis. It should be stressed that p-values are often misinterpreted in ways that lead to overstating the evidence against the null hypothesis when conventional thresholds (e.g., p-value = 0.05) that signify "statistical significance" are utilized [57]. As shown in [56], conventional levels of significance do not actually provide strong evidence against the null hypothesis. When the prior probabilities of the null and the alternative hypotheses are equal, the upper bound on the posterior probability of the alternative hypothesis equals 0.89 for a p-value of 0.01, which is often considered "highly significant"; hence, there remains at least an 11% chance that the null hypothesis is true. For a p-value of 0.005, the upper bound on the posterior probability of the alternative hypothesis equals 0.933, which reduces the chance that the null is true by about 50% relative to when the p-value equals 0.01 [56].

Mesh Convergence
The well-established Grid Convergence Index (GCI) method [58] was adopted in order to access the discretization error of all simulations presented herein. To this end, three simulations with constant refinement ratio r = 2 were performed for the 31A case of the Casson model with Parabolic profile. Denoting f 1 , f 2 , f 3 the solutions for the fine, medium, and coarse mesh respectively of a parameter of interest, the order of convergence p equals [58] The GCI between two successive mesh refinements i and j can then be formulated as where the safety factor F s is usually set equal to 1.25 [59], and e ij is the relative error of f. Finally, convergence is achieved when the ratio CGI 23 /(r p · CGI 12 ) satisfies the following condition CGI 23 r p · CGI 12 → 1.
Convergence was accessed by examining the three hemodynamic variables of interest, namely TAWSS, OSI, and RRT. The final mesh consisted of 776,000 first-order hexahedral elements (reference mesh). All simulations were performed using meshes with a total number of elements that do not vary more than 5% with respect to the reference mesh. Table 3 presents the corresponding results; one can observe convergence rates well within the acceptable regime for grid convergence.  Table 4 summarizes the findings of 14 robust, trimmed-mean-based ANOVA models, which evaluate the effects of alternative inlet velocity distributions and rheological models regarding the hemodynamic behavior of real aneurysm geometries. It is worth highlighting that the observed p-values clearly suggest that interaction effects are negligible; hence, the effects of IVD do not depend on the selected rheological specification and vice versa, for all response variables examined. Starting with flow asymmetry (Table 5), there is no evidence of different average levels across IVD or R m groups. Figure 5a supports this finding, with substantial overlapping in the group-specific boxplots. Table 4. p-values derived from heteroscedasticity-and outlier-robust, trimmed-mean-based ANOVA analyses for the effects of alternative inlet velocity distributions and rheological models at the hemodynamic behavior of real aneurysm geometries. Response variables related to hemodynamic behavior are shown in rows with the two main factors and their interaction in columns. p-values very close to zero, shown in bold, provide strong evidence against the null hypothesis of equal means across factor groups.   On the other hand, observed flow dispersion percentages differ across IVD specifications (Figure 5b), with the corresponding post hoc comparisons suggesting that it is the Parabolic inlet velocity specification that produces lower average levels for FDpct. Specifically, the difference in FDpct trimmed-means equals −6.982 for Parabolic vs. Plug (95% CI [−12.674, −1.289]); the corresponding p-value equals 0.012, which provides weak evidence against the null hypothesis of equal average levels. A noteworthy, although also weak, On the other hand, observed flow dispersion percentages differ across IVD specifications (Figure 5b), with the corresponding post hoc comparisons suggesting that it is the Parabolic inlet velocity specification that produces lower average levels for FD pct . Specifically, the difference in FD pct trimmed-means equals −6.982 for Parabolic vs. Plug (95% CI [−12.674, −1.289]); the corresponding p-value equals 0.012, which provides weak evidence against the null hypothesis of equal average levels. A noteworthy, although also weak, dissimilarity is also observed between Parabolic and Womersley specifications: the difference in FD pct trimmed-means equals −8.314 (95% CI [−15.626, −1.003]) and the corresponding p-value = 0.015.

ANOVA Models
Regarding thrombus-prone percentage areas (Table 6), it is interesting to observe that TAWSS-based assessments depend weakly on the selected rheological model, whereas OSI-based assessments depend strongly on the chosen inlet velocity distribution ( Figure 6). The post hoc comparisons suggest that the observed differences for TAWSSpct are mainly due to N vs. P and HB vs. P (N vs. P difference in TAWSSpct equals 14

Discussion
The above results carry several implications for methodologies developed to answer clinical questions. Computational modeling should provide insights regarding the physiology of the circulatory system, pathophysiology of cardiovascular diseases, and performance of vascular therapies. If one aims to predict thrombus deposition in order to characterize the risk profile of a given AAA, such a prediction needs to be accurate in order to be relevant from a clinical perspective. In other words, if therapeutic decisions and management of these lesions are going to be determined based on computational modeling and numerical simulations, these must provide sound rupture risk estimation.
Hemodynamic variables may exert direct clinical implications during AAAs evaluation, related to thrombus-prone regions and relevant tangential forces. Indeed, since intraluminal thrombus has been previously shown to be a significant determinant of the risk of rupture, the choice of the inlet profile or rheological model may affect AAA assessment. An investigation of the literature suggests that different rheological models and inlet velocity distribution profiles are used interchangeably in computational modeling, with uncertain implications for the results obtained. Deviations between different model assumptions should be recorded, and reporting standards regarding the methodology of such calculations should be developed.
The null hypotheses examined here essentially imply that practitioners with different choices regarding their simulation setups ought to produce vascular flows that do not differ substantially in terms of their main characteristics: TAWSS, OSI, percentages of thrombusprone regions, flow asymmetry, and dispersion. The outlier-and heteroscedasticity-robust statistical tests presented in Section 3 (which are conservative, relative to conventional, non-robust procedures) suggest that RRT, a variable of primary interest as it takes into account both the average tangential force on the wall over the entire cycle and the oscillatory nature of the flow, does not differ significantly across rheological models and inlet velocity distribution profiles. The same finding is also observed when one focuses on flow asymmetry, whereas differences in TAWSS levels are mainly due to the Newtonian specification; non-Newtonian alternatives that are widely applied in practice do not lead to strong evidence against the null hypothesis of equal levels.
Interestingly, our analyses suggest that the selection of an inlet velocity profile appears as a more critical factor relative to the rheological specification in vascular flow simulations. Furthermore, the two factors exert additive effects on variables that characterize vascular flows, with negligible interaction effects. Indeed, flow dispersion averages derived using the Parabolic profile are substantially lower relative to the ones corresponding to Plug or Womersley specifications. On the other hand, observed differences are statistically indistinguishable for alternative rheological models. If one focuses on thrombus-prone regions and the oscillatory characteristics of the flow (via OSI), it is mainly the Plug specification that produces substantially lower percentages relative to Parabolic and Womersley profiles.
Taking into account previous studies, one could assume that a Parabolic or Womersley inlet profile approximates patient-specific flow conditions more closely compared to the Plug profile, at least for critical hemodynamic variables such as thrombus-prone AAA segments [18,24]. Conversely, focusing on other variables such as flow dispersion, other patterns can be identified, such as the Parabolic profile being the specification that produced heterogeneity. Nevertheless, the relevance of this metric and the clinical implications of large versus small flow dispersion are not yet well understood. Patient-specific inflow velocity profiles might be advantageous, but these would require advanced imaging modalities and analysis, which can be challenging, impractical, and not always available.
As with the majority of studies, the present work is subject to limitations. The first one is the rigid wall assumption, neglecting the effect of wall compliance. A second limitation is the implementation of only one outlet pressure boundary condition. It should be noted though that due to the large number of rheology models and patient-specific geometries, it was practically not possible to incorporate the above two additional factors since the required simulations would have been in the order of thousands. To overcome this obstacle, a future work will include fewer geometries and blood flow models, which according to the results of this study are less critical with respect to the choice of inlet velocity profile.

Conclusions
This manuscript evaluates the effects of alternative modeling choices, namely alternative rheological specifications and inlet velocity distributions, on simulated blood flows. The finite volume method was employed, focusing on seven patient-specific abdominal aortic aneurysms. Fourteen variables that characterize and summarize vascular flows were computed; for each case, five rheological specifications and three inlet velocity boundary conditions were utilized. The significance of observed differences, which result from different modeling choices, were evaluated to our knowledge for the first time, with modern and robust statistical procedures.
Our findings can be summarized as follows. The two examined factors, namely rheological models and inlet velocity distributions, exert additive effects; to put it otherwise, their interaction effects were negligible on all variables that characterize vascular flows. The selected inlet velocity profile is a stronger factor, relative to the chosen rheology, for variables that are associated with the oscillatory characteristics of blood flow. Interestingly, response variables that relate to the average tangential force on the wall over the entire cycle do not differ significantly across alternative factor levels, as long as one focuses on non-Newtonian specifications.
Author Contributions: Conceptualization, K.T., Y.K. and N.K.; methodology, K.T., Y.K., N.K. and C.V.I.; software, K.T.; writing-original draft, K.T., Y.K. and N.K.; writing-review and editing, K.T., Y.K., N.K. and C.V.I.; supervision, K.T. All authors have read and agreed to the published version of the manuscript. Informed Consent Statement: Informed consent was obtained from all subjects involved in the study.
Data Availability Statement: Statistical investigations are fully reproducible as all produced and analyzed data are available upon request; replications of outlier-and heteroscedasticity-robust statistical outcomes can be performed with R package WRS2.

Acknowledgments:
The authors would like to thank Antonios Karasavvidis, CFD Applications, Customer Service, BETA CAE Systems SA for many useful discussions on mesh building.

Conflicts of Interest:
The authors declare no conflict of interest.